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Abstract. Uncertainty quantification in aerodynamic simulations calls for efficient numerical methods to reduce 
computational cost, especially for the uncertainties caused by random geometry variations which 
involve a large number of variables. This paper compares five methods, including quasi-Monte 
Carlo quadrature, polynomial chaos with coefficients determined by sparse quadrature and gradient- 
enhanced version of kriging, radial basis functions and point collocation polynomial chaos, in their 
efficiency in estimating statistics of aerodynamic performance upon random perturbation to the 
airfoil geometry which is parameterized by independent Gaussian variables. The results show that 
gradient-enhanced surrogate methods achieve better accuracy than direct integration methods with 
the same computational cost. 

Key words, aerodynamics simulation, airfoil geometric uncertainty, surrogate modeling, gradient-enhanced 
kriging, numerical integration 

AMS subject classifications. 65D05, 65D30, 65D32 

1. Introduction. Aerodynamic simulations are subject to uncertainties in aircraft geome¬ 
tries due to various unpredictable perturbations, e.g. manufacturing tolerances, operational 
wear and surface icing etc [11]. These uncertainties have the potential to dramatically lower 
aerodynamic performance. Therefore enabling uncertainty quantification is crucial for the 
robustness of aircraft designs. 

Uncertainty quantification (UQ) is a vivid, fast growing research field in the CFD com¬ 
munity. In the context of aerodynamic design, most UQ papers focus on the treatment of 
uncertainties in flow conditions, e.g. Mach number or angle of attack, see [36, 3, 40, 22], and 
only very little attention has been paid to uncertainties in the geometry. 

The quantification of geometry-induced uncertainties is computationally challenging due 
to large number of variables so that it calls for efficient numerical methods. Methods sample 
on a regular grid are usually less favoured due to two reasons. The first is their lack of 
tolerance to sample failures which are not rare in CFD computations. The second is their 
vulnerability to “curse of dimensionality”, though sparse grid techniques like those employed 
in [39, 25, 37] can alleviate the curse to some extent. Probabilistic collocation method based 
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on Gaussian quadrature is used in [30] and multivariate Gaussian quadrature in [8] to quantity 
airfoil geometric uncertainties, both with a small number of variables due to the curse. On 
the other side, methods based on scattered sampling schemes provide the robustness against 
sample failures and more freedom in choosing sample number N. 

Monte-Carlo (MC) and quasi-Monte Carlo (QMC) quadrature are direct integration meth¬ 
ods using scattered sample scheme. The MC quadrature and its variance-reduced kin (e.g. 
Latin Hypercube) have dimension-independent error convergence rate 0(N~ 2 ), and QMC has 
the worst case rate 0(log d (N)N~ 1 ) with d the number of variables. Due to higher degree of 
sample uniformity the latter can be more accurate even with a large d. 

Surrogate-based methods are gaining more attention recently, e.g. radial basis functions or 
kriging in [28, 4, 31, 14] and polynomial chaos expansion in [21, 10, 24, 26, 37]. A non-intrusive 
Galerkin way of constructing surrogates is proposed in [13]. [14] shows that a kriging surrogate 
is better than plain Monte Carlo and Latin Hypercube methods in estimating mean value of 
a bivariate Rosenbrock function. In a recent paper [37] the authors made a comparison of 
two classes of sample-based polynomial surrogates, generalized polynomial chaos in its sparse 
pseudospectral form and sparse grid stochastic collocation method, in a test case of uncertain 
airfoil geometry. The authors also proposed an improved version of the latter method with 
dimension adaptivity. 

Surrogate-based methods seem to be of advantage when the gradients can be sampled 
at a relatively lower cost, because direct integration methods like MC and QMC cannot 
effectively utilize the gradients (augmenting samples generated by Taylor’s expansion are not 
statistically independent hence bring little benefit to the accuracy of the integration). In 
[21] a point-collocation polynomial chaos surrogate is built to facilitate shape optimization, 
reinforced by gradient information acquired by local sensitivity analysis which costs at most 
20% of that for response quantities. In case of CFD model, even cheaper gradients can be 
obtained by using an adjoint solver [5]. By the adjoint solver all the partial gradients can be 
sampled at the cost of about one CFD model evaluation, i.e. each gradient costs only 1/d of 
that for response quantities. 

However in searching the literatures on quantifying geometry-induced uncertainties we 
have the following impressions: 

1. In this field only one efficiency comparison of various UQ methods (not only versus 
Monte-Carlo method) is found in the forementioned paper [37] where the authors focus in 
surrogate methods based on sparse grid stochastic collocation and generalized polynomial 
chaos. 

2. Few papers address this problem with gradient-employing surrogate method which we 
consider the most efficient. 

3. No effort is made to base a UQ method comparison on “accurate enough” reference 
statistics. 

In this paper, we compare the efficiency of various UQ methods in quantifying performance 
uncertainties caused by geometric uncertainties of airfoil in two test cases. In the first test 
case the comparison includes five methods consisting of one direct integration, namely quasi- 
Monte Carlo (QMC) quadrature, and four surrogate-based methods: polynomial chaos (PC) 
with the coefficients computed by a sparse Gauss-Hermite (SGH) quadrature (as introduced 



Quantification of geometry-induced aerodynamic uncertainties 


3 


in [2]), gradient-enhanced kriging (GEK), gradient-enhanced radial basis functions (GERBF) 
and gradient-enhanced polynomial chaos (GEPC) method (see method descriptions in Section 
2). To make the comparison more reliable, the reference statistics (against which the error 
of estimated statistics is judged) is obtained by an integration of four million CFD samples 
and its accuracy justified sufficient by using multipartition method. In the second test case 
three methods are compared, namely QMC, plain kriging and GEK. The number of CFD 
samples utilized by the UQ methods is kept small (< 500) in the both test cases to make the 
comparisons more meaningful to industrial applications. 

For the geometric uncertainties in the test cases, we focus on smooth geometry deviations 
from the original shape as those typically arise during the manufacturing process. Other 
uncertainties or errors, such as those due to physical, mathematical and numerical modeling 
approximations, are not considered in this work. 

The plan for the rest of the paper is as follows: Section 2 briefly introduces the methods 
in the comparison. Section 3 sets up a test case based on an Euler model of an airfoil with 
uncertain geometry, and gives results of the comparison as well as a discussion. Section 4 does 
the same for a Navier-Stoke test case. Section 5 gives a summary of the paper. 

2. Methods. In this section the methods in the comparison are introduced. More weights 
are given to two novel implementations, the gradient-enhanced radial basis functions and 
gradient-enhanced polynomial chaos method. Only brief introductions or references are given 
to the others. 

For the Quasi-Monte Carlo (QMC) quadrature [7], Sobol low discrepancy sequence gener¬ 
ated by the algorithm given in [19] is used since the sequence is shown comparable or slightly 
better than other sequences in a numerical comparison [35]. 

The gradient-enhanced kriging (GEK) and plain kriging [9] are implemented by using 
Surrogate-Modeling for Aero-Data Toolbox (SMART) [15] developed at DLR, opting for ordi¬ 
nary kriging and a correlation model of cubic spline type which is considered the most efficient 
in similar situations in [27]. The internal parameters of the correlation model are fine-tuned 
to fit the sampled data by a maximum likelihood estimation [45]. 

The gradient-enhanced radial basis function method is introduced in Section 2.1, poly¬ 
nomial chaos method based on a sparse Gauss-Hermite quadrature (PC-SGH) and gradient- 
enhanced polynomial chaos method (GEPC) are introduced in Section 2.2. 

To quantify uncertainty via surrogate methods, one first establishes the surrogate based 
on a small number of QMC samples of the CFD model, and integrates for the target statistics 
and probability density function (pdf) by a large number (we take 1 x 10 6 ) of QMC samples 
on the surrogate, only except for the polynomial chaos methods where the mean and variance 
can be directly obtained. 

2.1. Gradient-enhanced radial basis functions (GERBF). Denoting the CFD model as 
/, a radial basis function (RBF) [6] approximate takes the form 

N 

/(o = 

2=1 
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where fa are radial basis functions, || • || denotes the Euclidean norm, and is a sample point 
where fa is radial about. The coefficients Wi are determined by fitting /(£) to N samples. 

Denoting the Euclidean distance from the center as r, popular types of far) include 
y/r 2 + o 2 (multiquadric), \/y/r 2 + a 2 (inverse nrultiquadric), exp(—a 2 r 2 ) (Gaussian) and r 2 ln(ar)| 
(thin plate spline), in which a is a parameter to be fine-tuned for a particular set of samples. 
Gradient-assisted RBF were proposed in [12, 33] where the gradients are exploited in the same 
way as gradient-enhanced kriging (GEK). This method is numerically equivalent to a GEK 
that uses RBF as the correlation function. 

A different gradient-employing RBF method is implemented in this work, which only 
requires RBF to be first-order differentiable, in contrast to the second-order differentiability 
required by the gradient-assisted RBF. We call it gradient-enhanced RBF (GERBF) in this 
paper. To accommodate the gradient information, this method introduces additional RBF 
that are centred at non-sampled points, i.e. an GERBF approximation is 

N{l+d) 

/(«)= £ Wi fa(\\£ — £w||), with d the number of variables 
i =1 

The with i < N are sampled points, while those with i > N are non-sampled points which 
can be chosen arbitrarily as long as none of them duplicates (or too close to) the sampled 
ones. The coefficients w = {mi, • • • , ic’jv(n-rf) are determined by fitting /(£) to N samples 
and Nd gradients, i.e. solving the following system, 

Am = f. 

Abbreviating fa( ||£ — £^||) as fa(£), A is a iV(l + d)-sized square matrix whose first N rows 
contains values of fa(£) evaluated at the N sampled points and the rest rows values of dfa/d^j 
for 1 < j < A T evaluated at the N sampled points, f is a N(1 + d)-length vector containing 
the sampled values of / and its d partial derivatives. 

Since A is often ill-conditioned the system is solved through a truncated singular value de¬ 
composition [34, Chapter 15] in which singular values smaller than a threshold are discarded 1 . 
This gives a solution with the smallest L 2 norm which is less spoiled by the ill-conditionedness 
[34, Chapter 15]. For large systems a more efficient and stable solution by pivoted Cholesky 
decomposition is proposed in [29]. 

Inverse multiquadric RBF is used in this comparison. The internal parameter a is fine- 
tuned by a leave-one-out cross-validation as in [4]. 

2.2. Polynomial chaos methods. According to Wiener [43], a surrogate could also be 
made in the form of a truncated polynomial chaos expansion (PCE) : 

K 

(2.1) /(o = E c ^) 

i=0 

where H. L is a multivariate Hermite polynomial with Gaussian variables (see introduction in 
e.g. [32]). It is called “truncated” because the K is taken finite for practical computation. 

1 Let A i be the singular values, those |Ai|/maxi(|Ai|) < N( 1 4- d) x 10~ 13 are discarded. 
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K is related to the maximal order of polynomial (p) and number of variables (d) by K = 
(p + d)\/(jp\d\). 

The coefficients c* in Equation (2.1) could be determined by various ways[23, 44], in this 
work the following two of them are implemented. 

2.2.1. Projection based on sparse Gauss-Hermite quadrature (SGH). Due to the mu¬ 
tual orthogonality of Hi s with respect to their inner product, the coefficients c* can be ex¬ 
pressed by 

1 fH?(t)p(t) d£ 5 

where p(£) is Gaussian probability density. The numerator is approximated by using sparse 
Gauss-Hermite (SGH) quadrature, the denominator is computed analytically [32]. 

2.2.2. Gradient-enhanced point-collocation polynomial chaos method (GEPC). Point- 
collocation polynomial chaos method [18] determines q by fitting equation (2.1) to N collo¬ 
cation points. The work [17] suggests to use N = 2 K for the best performance which leads 
to an over-determined system to be solved by a least square approach. A gradient-employing 
version of this method is implemented in this work, in the name gradient-enhanced polyno¬ 
mial chaos (GEPC) method, where accordingly one chooses N{1 + d) = 2 K. The unknown 
coefficient c = {cq, c\, ■ ■ ■ ,ck} t is determined by solving the following system 


*c = f 

H 0 ^ W ) ••• H k (£&)- 

H 0 (^ N) ) ••• H k (^ n )) 

■■■ H%\^)' 

... H%>(£<"))_ 

with H^ = dHi/d^j. Thus, T' consists of values of the Hi and dHi/d^j evaluated at N 
sample points, sized N(1 + d) x K. The vector f is composed as in the GERBF method. 

Once a surrogate based on polynomial chaos is established, the estimates of the mean and 
the variance of /(£) can be directly obtained from c: 


in which 

= 


<S> 0 


= 


vR, = 

d>j> 0 


H ( 0 j \^) 


£ o 

rO) 




( 2 . 2 ) 


K r r 

p = c 0 , a 2 ° 2 i H h£)p{£)&£. 

i=l L ^ 


where the integration can be made analytically. The exceedance probabilities and pdf are 
integrated by a large number (we take 1 x 10 6 ) of QMC samples on the surrogate model. 
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3. An Euler test case of uncertain airfoil geometry. The first test case is a CFD model 
of inviscid flow around a RAE2822 airfoil at Mach number 0.73, angle of attack (AoA) 2.0 
degree and temperature 273.15 degree. We use flow solver TAU [16], opting for central flux 
discretization, scalar dissipation, backward Euler solver and “4w” multigrid cycle. The domain 
is discretized by a 193 x 33 grid in which the airfoil has 128 surface nodes, as shown in Figure 1. 



Figure 1. Grid for the RAE2822 airfoil in the Euler test case: the total grid (left) and zoom around the 
airfoil (right). 


The source of uncertainty is a random perturbation to the original airfoil geometry, which 
is modeled by a random field. Representing the original geometry (see Figure 4) by a matrix 
G £ M 128x2 w ] 1 j c } 1 contains the x,y coordinates of the surface nodes, the random field of 
perturbation can be written as R(g,ix) € M 128 with g = {x,y} € G and u € fl, 11 is the 
sample space of the random perturbation. 

R(g : u) is defined as a transformation of a Gaussian random field ip(g,uj) which can be 
easily parameterized by independent random variables by using Karhunen-Loeve expansions. 

The Gaussian field ip(g,u) is defined by a zero mean and a Gaussian-type correlation 
function: 

(3-1) co v(gi,gj) = a(gi) a(g :j ) exp ^ 9t ^ 


with correlation length l = 0.005, the standard deviation a is assumed to be: 


°( 9 ) 


(0.8-x) 0 - 75 if x < 0.8 
0 otherwise 


The zero perturbation to the trailing end is to guarantee the convergence of the CFD solver. 

To ensure boundedness of the perturbation, the random perturbation field R is made to 
have a sine-shaped probability distribution by such a transformation of the Gaussian field i/> : 


(3.2) 


R(g , u) = 2 s (. t ) • arccos (1 — 24> (ij) (g, co))) /it — l 
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where $ is the Gaussian cumulative distribution function and s(x) = (0.8 — x) 0 75 • \/0.00002. 

Perturbed geometries are generated by imposing the perturbation R to the original ge¬ 
ometry, in the direction normal to the airfoil surface: 

G(w) = G + R(g, uj) ■ n(g), 

where n(g) is the normal vector at g. 

However the above representation has 128 correlated variables while one prefers a repre¬ 
sentation with independent variables in a fewer number. A model reduction can be made by 
an approximation of if through a Karhunen-Loeve expansion (KLE) [1], i.e. 

k 

(3.3) i/>(g,u) ~ Vi{g) &(w) =: rjt(g,u) 

i—1 

where \ and V * denote the eigenvalues and eigenvectors of the covariance matrix generated by 
Equation (3.1) and G’s are independent standard Gaussian random variables. For our problem 
the number of variables is reduced to k = 9, this truncated KLE approximation retains 99.98% 
of variance of the random field. The KLE is optimal in all linear form representations with 
the same number of terms since it minimizes the loss in the variance caused by the model 
reductionfl]. In numerical implementation ijj(g,u>) is used instead of if)(g,uj) in Equation 
(3.2), so that the random field of perturbation R is parameterized by 9 independent Gaussian 
variables. In [20, 29] efficient algorithms for eigen-decomposition of large covariance matrices 
are given, by which one has the choice to compute only the k largest eigen modes. 

Figure 2 and 3 display three realizations of the i/> and R respectively. Examples of per¬ 
turbed geometries G are shown in Figure 4. Figure 5 shows the fast decline of the eigenvalues 

A*. 


Upper part Lower part 




Figure 2. Realizations of the Gaussian random field if: upper (left) and lower (right) part of airfoil. 

The test case is a scaled-down model with an inviscid flow assumption and relatively 
coarse mesh. This simplifying configuration is to facilitate a large number of evaluations 
to obtain reliable reference statistics of performance on the basis of which the accuracy of 
estimated statistics are compared. Since this accuracy is generally independent of the model 
bias caused by the particular configuration, the result of this comparison holds meaningful and 
representative in judging the efficiency rank of UQ methods in other stricter configurations. 

3.1. Settings of the comparison. With the above parameterization of geometric uncer¬ 
tainties, the five UQ methods are applied to the test case and compare their efficiency in 
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Figure 3. Realizations of the transformed random field R: upper (left) and lower (right) part of airfoil. 



Figure 4. Three realizations of perturbed geometry 


Decline of eigenvalues 



modes 


Figure 5. Decline of eigenvalues 
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estimating statistics of two aerodynamic system response quantities (SRQ), namely the co¬ 
efficients of lift and drag (Cl and Cd), as well as their probability distribution functions 
(pdf). 

3.1.1. Target statistics and their reference values. The following statistics of Cl and 
Cd are estimated: 

• means hl and ft d 

• standard deviations er l and a q 

• exceedance probabilities Pl )K = Pro {Cl < ftp ~ k ■ &l} and Pd,k = Pro {Cd > 
Hd + k • &d} with k = 2,3 

The reference values of these statistics are obtained from a relatively large number (N = 
4 x 10 6 ) of quasi-Monte Carlo (QMC) samples of the CFD model. To justify the validity of 
using these reference values in the efficiency comparison, their accuracy are estimated using 
Snyder’s multipartition method [41] since the theoretical error bound of QMC integration 
is not a practical accuracy indicator. Suppose S denotes any of the above statistics, by 
this method one makes an equal-size ?n-partition of the 4 x 10 6 samples, and obtained m 
estimates by integrating on each partition only, then compute the sample standard 

deviation of S t . is an estimate of the standard deviation of QMC estimate of S with 
sample number N/m. To extrapolate to ci one computes for four values of m, i.e. 
m = {128,64,32,16}, and fit a line across the four (\og(N/m), log(c m )) points using weighted 
linear least square method, the m values are used as weights to account for the increasing 
variability for smaller values of m. With the fitted slope a and intercept (3 the ci is extrapolated 
as ci = exp(a log(lV) ±/3) which is an estimate of the standard deviation of the QMC estimate 
of S with N samples. 

values for each reference statistics are listed in Table 1. In the efficiency comparison 
of the UQ methods (results given in section 3.2) the smallest measured errors (measured 
against these reference values) in mean and standard deviation (stdv) are at least by 10 
times larger than 3 x <q, which means, by taking the assumption that the reference values 
are Gaussian distributed around the true values of the statistics, the measured errors have a 
99.73% confidence interval of at widest ±10%. For exceedance probabilities this confidence 
interval is also valid except for a few measured errors, as shown in Figure 6 and 7 where the 
values of the corresponding 3xq are depicted by thick dash line. 


?i (ul) 

?i(ol) 

Cl(-Pi,2) 

Ci(Tl,3) 

9.1e-9 

5.0e-8 

1.8e-5 

9.3e-6 


?i(/uj) 

?i(o- D ) 

Ci (Pd,2) 

Ci {Pd,3) 

3.0e-9 

1.4e-8 

1.2e-5 

6.8e-6 


Table 1 


Estimated stdv of reference statistics for Cl (top) and Cd (bottom) 


3.1.2. Measure of cost. The computational cost is measured in term of “compensated 
evaluation number” N c . For the three gradient-employing methods, N c = 3 N with N denoting 
the number of CFD model evaluations, since with an adjoint solver the computational cost 
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of the gradients of each system response quantity (SRQ) equals approximately to the cost of 
one CFD evaluation and we need the gradients for two SRQs (Cl and Cd)- For QMC and 
PC-SGH methods, N c = N. 

3.1.3. Range of sample numbers . The number of CFD evaluations is kept smaller than 
500 in this comparison, as shown in the Table 2. Most of the methods share the same range 
except for PC-SGH in which only quadrature level 2 and 3 are used since the sample number 
of level 4 is well beyond 500. 


Methods 

N 

N c 

QMC 

30 ~ 420 

30 - 420 

PC-SGH (order 2, 3) 

19 and 181 

19 and 181 

GEK 

O 

rH 

l 

O 

30 ~ 420 

GERBF 

10 ~ 140 

30 - 420 

GEPC (order 2,3 and 4) 

11, 44 and 143 

33, 132 and 429 


Table 2 

Sample number in the comparison 


3.2. Numerical results and discussion . The results of the efficiency comparison are 
shown in Figures 6-9. Figures 6 and 7 show the errors of the five methods in estimating the 
target statistics of Cl and Cd- It is observed there that generally the gradient-employing 
surrogate methods perform better than direct integration methods. This can be partially 
ascribed to that the former utilize more information with the same computational cost N c , 
i.e. they use jt^N c conditions with s the number of SRQs, while a direct integration method 
uses N c conditions. This advantage comes from the cheaper cost of the gradients computed 
by an adjoint solver when s is smaller than d (in this test case, 2 versus 9), and the advantage 
would increase for a larger d. 

The PC-SGH method has only two data points due to the very limited choices of sample 
number. It is hard to evaluate its error convergence property on only two data points. This 
shows a shortcoming of sparse grid quadratures with a relative high number of variables. 
The GEPC method shows constant convergence in error, but not as efficient as the other 
two gradient-employing methods. The reason could be that the polynomial surrogate tend to 
“overshoot” in the outskirt of the domain. But the PC methods have a merit that they do 
not need a parameter optimization procedure. A future combination of gradient-employing 
and the more sophisticated PC methods introduced in [37] has the potential to enhance the 
efficiency. 

GEK and GERBF are the most efficient methods as far as seen from these results. Beside 
the cheaper gradients, this could be attributed to properties of the kernel functions they 
use and the effort of tuning parameters. Convergence rate of inverse multiquadric RBF was 
estimated 0 (e _7//ft ) with h the fill distance and 7 a constant, which translated to a rate 
in N at 0(e~' yNl/d ) [42]. We assume the error in the statistics is proportional to that in 
point-wise approximations, and find in this 9-variate test case the observed convergence rate 
is much better than 0(e~' yNl/9 ), this hints that the “effective” fill distance h reduces faster 
than 0{N~ l / d ) due to that some variables are less important than the others (typical for a 















Quantification of geometry-induced aerodynamic uncertainties 


11 


On estimating g L 


On estimating o L 



On estimating P L 



On estimating P L3 



Figure 6. Absolute error in estimating mean, standard deviation (upper row) and exceedance probabilities 
(lower row) of Cl 


KLE-based parameterization). This is exploited by anisotropic kernel functions out of the 
parameter tuning. 

GEK seems slightly better than GERBF according to the results, especially at smaller 
N c values. But this might not be generalized to other applications. The advantage of GEK 
could come from the possible advantage of cubic spline kernel over inverse multiquadric in this 
particular case. Different ways of utilizing gradients by GEK (involving 1st and 2nd-order 
derivatives and generating a symmetric kernel matrix) and GERBF (involving only 1st order 
derivative and generating non-symmetric matrix) could also contribute to the difference in 
their performance. 

At larger N c values a “rebound” of error can be observed in GERBF and GEK results 
which is caused by numerical instablities due to high condition number of their system ma¬ 
trices, in spite of the stabilizing treatments 2 set a priori in the two methods. Though there is 
an uncertain relation (or trade-off principle) exists stating that accuracy and stability cannot 


2 For the GERBF, the truncated SVD as introduced in section 2. For the GEK implemented by SMART 
the system matrix is regularized by adding a 10 -13 “nugget” to the diagonal. They could be insufficient when 
size of the system matrices become large enough. 
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On estimating n D On estimating o D 





Figure 7. Absolute errors in estimating mean, standard deviation (upper row) and exceedance probabilities 
(lower row) of Cd 


be good at the same time [38], more advanced stabilizing techniques, e.g., pivoted Cholesky 
decomposition [29] could improve the convergence as sample number is large. 

For the mean and standard deviation the reliability of the measured errors are justified by 
that they are at least by 10 times larger than the corresponding 3 x <q (<q defined in 3.1.1). 
For exceedance probabilities the 3 x cl values are not that small (as indicates by the thick 
dash lines in the pictures) but are still small enough for most of the measured errors. 

Figures 8 and 9 show the probability density functions (pdf) of Cl and Cd estimated 
by QMC and the three gradient-employing surrogate methods with N c = 33 (the smallest 
possible N c value for GEPC method), comparing with the reference pdf (computed through 
4x 10 6 QMC samples). There one observes that for the same computational cost, the surrogate 
methods yield much more accurate pdf’s. This is consistent with their relative performance 
in estimating the statistics. 

However, the advantage of the gradient-employing surrogate methods should not be taken 
universal. For example, on problems with few variables or that are mostly undifferentiable the 
advantage might not exist. If the problem has a very oscillatory topography the advantage 
could only be manifested with relatively more samples. 
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pdf of C L by QMC, N=33 pdf of C L by GEK, N,=33 




100 


-reference pdf 

100 


-reference pdf 



-pdf by GEPC 



-pdf by GERBF 

80 

l/\ 

// \ 
i \ 


80 

lf\ 

l \ 

// \ 


60 

/ \ 


60 

// \ 
i \\ 

// \\ 


Q_ 

/ u 

/ \\ 


a. 

f \\ 

/ n 


40 

/ A 

/ A 

/ A 


40 

/ u 

/ \\ 

/ A 


20 

/ A 

/ n 

/ 


20 

/ A 

/ A 

/ \ 

/ 



y v 




\\ 

\\ 

n 

_ 

Vs- 

n 

— — 



o 1 -—— '-‘-‘- *--- o 1 —-■-■--■— 

0.79 0.8 0.81 0.82 0.83 0.79 0.8 0.81 0.82 0.83 

C L C L 


Figure 8. Estimated pdf (in dash line) of Cl by QMC, GEK, GEPC and GERBF at N c = 33. 

4. A Navier-Stoke test case of uncertain airfoil geometry. The second test case is also 
based on a RAE2822 airfoil, but here we use a Reynolds averaged Navier-Stokes(RANS) solver 
in the TAU package [16], opting for a turbulence model SAO, with a nominal Mach number 
0.729 and angle of attack(AoA) 2.31 degree. The domain is discretized by a much finer grid 
in which the airfoil has 380 surface nodes, as shown in Figure 10. 

The sources of uncertainty are random perturbation in the airfoil geometry and operational 
parameters Mach number and AoA. The geometry perturbation is modeled by a zero-mean 
random field through a Karhunen-Loeve expansions (KLE) using the same methodology as 
introduced in the first test case, but with a different setting. Here the upper and lower 
surfaces of the airfoil are treated as two separate random fields as their correlation in geometric 
variations is assumed weak in this test case. The standard deviation a in (3.1) is uniformly set 
to 1 in the fields. After the standard Gaussian random fields are approximated through the 
KLE in (3.3), they are transformed to perturbation fields R being of Beta distribution 3 within 
support [—6(g), 5(g)] for the boundedness. 5(g) is set to 1.5% of the airfoil thickness at the 

3 The probability density function is P(x) = ^ Xma - x ~ x '> ( ' <3 -- with [i m ™,i raal ] support of the 

distribution, B Beta function and a = /3 = 4. 
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Figure 9. Estimated pdf (in dash line) of Cd by QMC, GEK, GEPC and GERBF at N c = 33. 



Figure 10. Grid for the RAE2822 airfoil in the RANS test case : the total grid (left) and zoom around the 
airfoil (right). 
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node g. We keep all the eigenvalues larger than ICE 7 in the truncated KLE approximations, 
this leads to a parametrisation in 24 variables for the geometric perturbation only. Figure 11 
display four realizations of the and R respectively. Examples of perturbed geometries G 
are shown in Figure 12. 






Figure 11. Four realizations of the Gaussian random field ip (left) and transformed field R (right). 

The perturbations in Mach number and AoA are of the same Beta distribution, but with 
a support within ±2% of the nominal values, adding up to a total number of variables as 26 
for the test case. 

Three UQ methods, quasi-Monte Carlo quadrature, gradient-enhanced kriging and plain 
kriging, are applied to the test case and their efficiency compared in estimating four statistics 
(mean, standard deviation, skewness and kurtosis) of the coefficient of lift (Cl). 

The accuracy of the estimates is judged by comparing to reference statistics which are 
based on 10000 QMC samples. 

The computational cost is again measured in terms of “compensated evaluation number” 
N c as introduced in the first test case. The only difference is that for the gradient-employing 
method (GEK) N c is set to 2 N since here we only handle one system response quantity (SRQ) 

C L . 

4.1. Numerical results and discussion . For this RANS test case our reference statistics 
at hands are not as accurate as those in the Euler one due to the much fewer CFD evaluations 
affordable. Therefore displaying error convergence in logarithmic scale is not fully supported 
by the precision. So we are content with observing the approaching of the estimated statistics 
to the references in a real scale. In Figures 13 we contrast the estimates of the statistics to 
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Figure 12. Four realizations of perturbed RAE2822 airfoil geometry, 10 times exaggerated. 

the reference statistics. As expected, generally the estimates approach to the references along 
the increase of the cost measure N c . It is observed that the gradient-enhanced kriging (GEK) 
outperforms plain kriging and quasi-Monte Carlo quadrature(QMC) in the estimation of all 
statistics. As in the first test case this can be explained by the more information utilized by 
the former method at the same computational cost. Especially the contrast of GEK and plain 
kriging highlights the efficiency gain brought by the cheaper gradient information. 

5. Conclusion. This paper compares efficiency of various methods in quantifying geometry-| 
induced aerodynamic uncertainties. An Euler and a RANS test cases based on RAE2822 
airfoil are set up by perturbing the geometry by random fields which are approximated and 
parametrized through Karhunen-Loeve expansions. UQ methods including Quasi-Monte Carlo 
quadrature, polynomial chaos with coefficients determined by sparse grid quadrature, gradient- 
enhanced radial basis functions, gradient-enhanced polynomial chaos, gradient-enhanced krig¬ 
ing and plain kriging are applied to the test cases and are compared in their efficiency in 
estimating some statistics and probability distribution of the uncertain lift and drag coeffi¬ 
cients. The results show that gradient-employing surrogate methods are more efficient. The 
advantage is due to the cheaper gradients obtained by using adjoint solver, and is expected 
to increase with an increasing d. 
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Figure 13. Estimates of mean, standard deviation, skewness and kurtosis of Cl, contrasted to the reference. 
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